In [17]:
conf_fig(-10)
Out[17]:
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>
In this notebook, we draw maps that measure, in some way, how solid are the predictions used for the horizontal sections in https://www.mdpi.com/2077-1312/10/7/986. For this purpose we take a grid of points in the Llobregat River Delta and we assign to each point a confidence metric satisfying some desirable properties.
First, we configure the directory where the data can be found and also the directory where the figures will be saved.
DATADIR='data/' # Directory with the data
FIGURESDIR='figures/' # Figures produced
Now, we import the packages and the functions that we are going to use. The customer functions are dfined in the file functions.py.
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import os
from sklearn.neighbors import NearestNeighbors
from sklearn.neighbors._base import _get_weights
import plotly.offline as go_offline
import plotly.graph_objects as go
import functions
from functions import *
The metric we propose uses the weight of the nearest neighbor and divies it by the sum of the weights of the four nearest neighbors such that each one belongs to a different class, this includes the nearest neighbor used to infer the granulometry, so the confidence degree lies always in the interval $[0.25,1]$.
A bonus of confidence is applied when there are several nearest neighbors of the class basement closer than any other non-basement neighbor.
# Global variables needed in the following function
# They are NearestNeighbors objects from sklearn that will be fit before calling the function
knn_sand = NearestNeighbors(n_neighbors=1)
knn_grav = NearestNeighbors(n_neighbors=1)
knn_clay = NearestNeighbors(n_neighbors=1)
knn_base = NearestNeighbors(n_neighbors=4)
The function confidence(X) will do the job.
def confidence(X): # The variable `X` is a row of points in the grid.
# For each class, we take the nearest neighbor of each point in X.
neight_grav=knn_grav.kneighbors(X)[0].reshape((len(X)))
neight_sand=knn_sand.kneighbors(X)[0].reshape((len(X)))
neight_clay=knn_clay.kneighbors(X)[0].reshape((len(X)))
# For the basement, we take up to 4 nearest neighbors to apply the confidence bonus
m = min(knn_base.n_samples_fit_, 4)
neight_base=knn_base.kneighbors(X, m)[0].reshape((len(X),m))
conf=np.zeros(len(X)) # Array that will store the confidence degree of each point
for i in range(len(X)): # Here, we conpute the confidence for each point
bas=sorted(list(neight_base[i]))
no_bas=[neight_grav[i],neight_sand[i],neight_clay[i]]
m_no_bas=min(no_bas)
w=_get_weights(np.array([no_bas+[bas[0]]]), 'distance')[0]
conf[i]=np.max(w)/np.sum(w)
# Bonus of confidence when several basement points are closer than the nearest non-basement neighbor.
s=1
while s<m and bas[s]<m_no_bas:
w=_get_weights(np.array([no_bas+[bas[s]]]), 'distance')[0]
conf[i]=conf[i]+(1-conf[i])*(w[3]-max(w[:3]))/np.sum(w)
s+=1
return 100*conf # We take the confidence to the range 25-100%
The function conf(data,poly,height,axi,grid), were the variables meaning are:
computes the confidence map at a given height. The output of this function is a list [xx,yy,C] where:
xx is the X UTM coordinate of a point in the grid.yyis the Y UTM coordinate of a point in the grid.C is a matrix where C[i,j] is the confidence estimated for the point (i,j).def conf(data,poly,height,axi,grid):
xx, yy = np.meshgrid(np.linspace(axi[0],axi[1],grid),np.linspace(axi[2],axi[3],grid)) # We create the grid
C=np.zeros(xx.shape,dtype = float) # We initalize the matrix
outside=[] # Points outside the contour
for i in range(xx.shape[0]):
for j in range(xx.shape[1]):
if poly.distance(geometry.Point(xx[i,j],yy[i,j]))>1:
outside.append((i,j))
df=data[data.Cota==height] # We take the data points at the given height
# We split the data into the four categories
gravels=df[df.Clase=='gravillas y gravas']
sands=df[df.Clase=='arenas']
clays=df[df.Clase=='arcillas y limos']
basements=df[df.Valor=='S']
# We fit the global NearestNeighbors objects
knn_grav.fit(list(zip(gravels['UTM_X'],gravels['UTM_Y'])))
knn_sand.fit(list(zip(sands['UTM_X'],sands['UTM_Y'])))
knn_clay.fit(list(zip(clays['UTM_X'],clays['UTM_Y'])))
knn_base.fit(list(zip(basements['UTM_X'],basements['UTM_Y'])))
# We call the confidence function
for i in range(xx.shape[0]):
d=list(zip(xx[i],yy[i]))
C[i]=confidence(d)
for (i,j) in outside: # We set the confidence to nan outside the contour
C[i,j]=np.nan
return [xx,yy,C]
We decide to create a confident map every 5 meters depth.
Heights=list(range(0,-101,-5)) # Heights where the confidence maps will be obtained
We also set the plot frame.
FRAME=500 # Plot frame
We read the data.
boreholes_data=pd.read_excel(DATADIR+'hsd new basements.xls') # We read the data from the xls sheet
boreholes_data=boreholes_data.drop(columns=['Codi']) # We do not need this column
ax=[boreholes_data.UTM_X.min()-FRAME,
boreholes_data.UTM_X.max()+FRAME,
boreholes_data.UTM_Y.min()-FRAME,
boreholes_data.UTM_Y.max()+FRAME] # Axis limits (with frame)
minx_rounded = 1000 * round(ax[0]/1000) # We round the X coordinate limits
maxx_rounded = 1000 * round(ax[1]/1000)
contour=pd.read_csv(DATADIR+'deltacontour.csv')
contour=contour.drop(columns=['Cota'])
contour_poly=geometry.Polygon(zip(contour['UTM_X'],contour['UTM_Y'])) # Delta polygon
xpol,ypol = contour_poly.exterior.xy # Delta contour coordinates
xx,yy,C=conf(boreholes_data,contour_poly,-100,ax,200) # Confidence map at 100m b.s.l. (below sea level)
The function coordinates creates a list with three lists (UTM_X, UTM_Y and Height coordinates).
xyzcontour=coordinates(DATADIR+'deltacontour.csv',[0,1,2])
bound=bounds(xyzcontour) # bounds for the figures
The funtion conf_fig(height) creates a figure with the confidence index at a given heigh.
def conf_fig(height):
# Confidence ranges
cmap = plt.cm.rainbow
cmaplist = [cmap(i) for i in range(cmap.N)]
cmap = mpl.colors.LinearSegmentedColormap.from_list('Custom cmap', cmaplist, cmap.N)
bounds=np.array([25, 30, 35, 40, 50, 60, 70, 80, 90, 95, 100])
norm1=mpl.colors.BoundaryNorm(bounds, cmap.N)
# figure data
xx,yy,C=conf(boreholes_data,contour_poly,height,ax,200)
# create de figure
fig, ax1 = plt.subplots()
ax1.imshow(C, cmap=cmap, origin='lower', norm=norm1, alpha=0.6,
extent=[xx.min(), xx.max(), yy.min(), yy.max()])
sm = plt.cm.ScalarMappable(norm=norm1, cmap=cmap)
sm.set_array([])
fig.colorbar(sm,spacing='proportional').set_label(' %', rotation=0)
ax1.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5)
ax1.set_xlabel('UTM_X')
ax1.set_ylabel('UTM_Y')
ax1.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))
ax1.set_title("Height "+ str(height) +' m')
ax1.axis(ax)
return ax1
For example, the confidence map at height -10 is
conf_fig(-10)
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>
The function knn_fig(h) uses the funtion knn we defined in the paper https://www.mdpi.com/2077-1312/10/7/986, which can be found explicitly in function.py, to draw a figure with the knn precictions at height h.
def knn_fig(h):
# figure settings and bounds: we set a color for each category
SCATTER=True
OPACITY=0.6
cm = mpl.colors.ListedColormap(['cyan','black','yellow','black','gray','black','saddlebrown'])
cm = cm(np.arange(cm.N))
cm[:,-1]=np.array([OPACITY,0,OPACITY,0,OPACITY,0,OPACITY])
cm = mpl.colors.ListedColormap(cm)
epsilon=10**(-15)
levels = [-0.5,0.5,1-epsilon,1+epsilon,1.5,2+epsilon,2.5,3.5]
norm=mpl.colors.BoundaryNorm(levels, 7)
# figure data
layer=layer_function(boreholes_data,h)
xyC=knn(boreholes_data,contour_poly,h,ax,200)
# create de figure
fig2, ax2 = plt.subplots()
ax2.contourf(xyC[0],xyC[1],xyC[2], levels, cmap=cm, norm=norm, antialiased = True)
if SCATTER:
ax2.scatter(layer[1]['UTM_X'], layer[1]['UTM_Y'], c='blue', s=5, alpha=0.8, label='gravels')
ax2.scatter(layer[2]['UTM_X'], layer[2]['UTM_Y'], c='orange', s=5, alpha=0.7, label='sands')
ax2.scatter(layer[3]['UTM_X'], layer[3]['UTM_Y'], c='black', s=5, alpha=0.6, label='clays')
ax2.scatter(layer[4]['UTM_X'], layer[4]['UTM_Y'], c='brown', s=5, alpha=0.6, label='basement')
ax2.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5)
ax2.set_xlabel('UTM_X')
ax2.set_ylabel('UTM_Y')
ax2.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))
ax2.set_title("Height "+str(h)+' m')
ax2.axis(ax)
return ax2
For example the knn prediction figure at height -10 is
knn_fig(-10)
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>
We define the funtion cyk_fig(h) to create a figure with the knn predictions and the confidence maps together at a given height h.
def cyk_fig(h):
cmap = plt.cm.rainbow
cmaplist = [cmap(i) for i in range(cmap.N)]
cmap = mpl.colors.LinearSegmentedColormap.from_list('Custom cmap', cmaplist, cmap.N)
bounds=np.array([25, 30, 35, 40, 50, 60, 70, 80, 90, 95, 100])
norm1=mpl.colors.BoundaryNorm(bounds, cmap.N)
SCATTER=True
OPACITY=0.6
cm = mpl.colors.ListedColormap(['cyan','black','yellow','black','gray','black','saddlebrown'])
cm = cm(np.arange(cm.N))
cm[:,-1]=np.array([OPACITY,0,OPACITY,0,OPACITY,0,OPACITY])
cm = mpl.colors.ListedColormap(cm)
epsilon=10**(-15)
levels = [-0.5,0.5,1-epsilon,1+epsilon,1.5,2+epsilon,2.5,3.5]
norm2=mpl.colors.BoundaryNorm(levels, 7)
# figure data
xx,yy,C=conf(boreholes_data,contour_poly,h,ax,300)
layer=layer_function(boreholes_data,h)
xyC=knn(boreholes_data,contour_poly,h,ax,300)
# figure creation
fig = plt.figure(figsize=(14, 4))
# Adds subplot on position 2
ax1 = fig.add_subplot(122)
# Adds subplot on position 1
ax2 = fig.add_subplot(121)
ax1.imshow(C, cmap=cmap, origin='lower', norm=norm1, alpha=0.6, aspect="auto",
extent=[xx.min(), xx.max(), yy.min(), yy.max()])
ax1.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5)
ax1.set_xlabel('UTM_X')
#ax1.set_ylabel('UTM_Y')
ax1.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))
ax1.set_title("Height "+ str(h) +' m')
ax1.axis(ax)
ax1.set_yticklabels([]) # Removes Y axis labels
ax2.contourf(xyC[0],xyC[1],xyC[2], levels, cmap=cm, norm=norm2, antialiased = True)
if SCATTER:
ax2.scatter(layer[1]['UTM_X'], layer[1]['UTM_Y'], c='blue', s=5, alpha=0.8, label='gravels')
ax2.scatter(layer[2]['UTM_X'], layer[2]['UTM_Y'], c='orange', s=5, alpha=0.7, label='sands')
ax2.scatter(layer[3]['UTM_X'], layer[3]['UTM_Y'], c='black', s=5, alpha=0.6, label='clays')
ax2.scatter(layer[4]['UTM_X'], layer[4]['UTM_Y'], c='brown', s=5, alpha=0.6, label='basement')
ax2.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5)
ax2.set_xlabel('UTM_X')
ax2.set_ylabel('UTM_Y')
ax2.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))
ax2.set_title("Height "+str(h)+' m')
ax2.axis(ax)
sm = plt.cm.ScalarMappable(norm=norm1, cmap=cmap)
sm.set_array([])
fig.colorbar(sm,spacing='proportional', ax=[ax2,ax1]).set_label(' %', rotation=0)
return (fig,ax1,ax2)
For example at height -10 we get
cyk_fig(-10)
(<Figure size 1008x288 with 3 Axes>,
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X'>,
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>)
We create a list with all figures every 5m.
figures=[cyk_fig(h) for h in Heights]
/var/folders/fd/bd2ntp3n1lq8rv9w3gj8d8gw0000gn/T/ipykernel_68207/739094612.py:26: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
And we save the figures in the figures directory.
for i in range(21):
filename='height'+str(-5*i)+'.png'
filename=FIGURESDIR+filename
os.makedirs(os.path.dirname(filename), exist_ok=True)
figures[i][0].savefig(filename, transparent=False, dpi=164, bbox_inches='tight')
The confidence for basement. We use a $200\times200$ grid.
# knn algorith predictions for a square grid of width 200
data_split=[[],[],[],[]]
for h in Heights:
cf=conf(boreholes_data,contour_poly,h,ax,200)
cc=knn(boreholes_data,contour_poly,h,ax,200)
for c in range(4):
x=[cc[0][i,j] for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
y=[cc[1][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
z=[h for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
ccf=[cf[2][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
data_split[c].append([x,y,z,ccf])
gravels_knn200=data_split[0] # knn predictions for gravels
sands_knn200=data_split[1] # knn predictions for sands
clays_knn200=data_split[2] # knn predictions for clays
basement_knn200=data_split[3] # knn predictions for basement
We group the confidence and the coordinates, getting four different lists.
x_bas200=np.concatenate([basement_knn200[i][0] for i in range(21)]) # X UTM coordinates
y_bas200=np.concatenate([basement_knn200[i][1] for i in range(21)]) # Y UTM coordinates
z_bas200=np.concatenate([basement_knn200[i][2] for i in range(21)]) # Z coordinates
conf_bas200=np.concatenate([basement_knn200[i][3] for i in range(21)]) # Confidence
xyzc_bas200=list(zip(x_bas200,y_bas200,z_bas200,conf_bas200)) # We zip the four lists
xyz_bas200=[x[:3] for x in xyzc_bas200] # Here, we take only the XYZ coordinates
We can work out the average confidence of the basement point for a map. Also the maximum and minimum heights where basement can be found and the medium height of the basement points.
bas_conf_mean=np.mean(conf_bas200)
bas_height_max=np.max(z_bas200)
bas_height_min=np.min(z_bas200)
bas_height_mean=np.mean(z_bas200)
(bas_conf_mean,bas_height_max,bas_height_min,bas_height_mean)
(77.98397866313435, 0, -100, -70.41003196387203)
Now, we set the grid size to 50, it is big enough and allows us to perform the computation within a reasonable time.
basementknn50max=pd.read_csv(DATADIR+'basementknn50max.csv')
basementknn50max
| UTM_X | UTM_Y | Cota | |
|---|---|---|---|
| 0 | 426262 | 4580059 | 0 |
| 1 | 426640 | 4580059 | 0 |
| 2 | 426640 | 4580345 | 0 |
| 3 | 427019 | 4580345 | 0 |
| 4 | 426640 | 4580632 | 0 |
| ... | ... | ... | ... |
| 11731 | 419445 | 4582063 | -35 |
| 11732 | 419824 | 4582063 | -35 |
| 11733 | 418688 | 4582350 | -35 |
| 11734 | 419067 | 4582350 | -35 |
| 11735 | 419445 | 4582350 | -35 |
11736 rows × 3 columns
basement50_coordinates=coordinates(DATADIR+'basementknn50max.csv',[0,1,2])
# We change the limits taking into account the data in basement_coordinates
basement50_bounds=[min(basement50_coordinates[0]),max(basement50_coordinates[0]),
min(basement50_coordinates[1]),max(basement50_coordinates[1]),
max(basement50_coordinates[0])-min(basement50_coordinates[0]),
max(basement50_coordinates[1])-min(basement50_coordinates[1])]
new_bounds50=bounds_join(bound,basement50_bounds)
We now interpolate
interpolation_knnbasement50=interpolation(basement50_coordinates,50,new_bounds50)
and reduce the data to the points inside the contour:
interpolation_knnbasement50cutted=cutting(interpolation_knnbasement50,contour_poly,500)
We shape the data as we need
xyz_interpolation_knnbasementcutted50=[[interpolation_knnbasement50cutted[1][i],
interpolation_knnbasement50cutted[2][i],
interpolation_knnbasement50cutted[0][i]] for i in range(len(interpolation_knnbasement50cutted[0]))]
We set the grid side to 200.
The results from the following four cells are stored in a csv file to save computation time.
Running them is not mandatory if the file basementknn200confidence.csv is present in the data folder.
The reader may just load the data as shown further.
# Data for the points where the basements is first reached
basementknn200max=pd.read_csv(DATADIR+'basementknn200max.csv')
# x,y,z list for the data in basementknn200max
basement200max_coordinates=coordinates(DATADIR+'basementknn200max.csv',[0,1,2])
# We change the limits taking into account the data in basement_coordinates
basement200_bounds=[min(basement200max_coordinates[0]),max(basement200max_coordinates[0]),
min(basement200max_coordinates[1]),max(basement200max_coordinates[1]),
max(basement200max_coordinates[0])-min(basement200max_coordinates[0]),
max(basement200max_coordinates[1])-min(basement200max_coordinates[1])]
new_bounds200=bounds_join(bound,basement200_bounds)
xyz_basement200max=list(zip(basement200max_coordinates[0],
basement200max_coordinates[1],
basement200max_coordinates[2]))
xyzc_basement200max=[x for x in xyzc_bas200 if ((int(x[:3][0]),int(x[:3][1]),int(x[:3][2])) in xyz_basement200max) ]
# We save the reduced data to a csv file to save computation time
df=pd.DataFrame(xyzc_basement200max,columns=["UTM_X","UTM_Y","Cota","Confidence"])
df.to_csv(DATADIR+'basementknn200confidence.csv',index=False)
Here, the csv file is loaded.
bknn200cmax=pd.read_csv(DATADIR+'basementknn200confidence.csv')
bknn200cmax
| UTM_X | UTM_Y | Cota | Confidence | |
|---|---|---|---|---|
| 0 | 426062.472362 | 4.579741e+06 | 0 | 28.593293 |
| 1 | 426155.718593 | 4.579741e+06 | 0 | 30.005940 |
| 2 | 425969.226131 | 4.579812e+06 | 0 | 28.314688 |
| 3 | 426062.472362 | 4.579812e+06 | 0 | 31.066904 |
| 4 | 426155.718593 | 4.579812e+06 | 0 | 33.898017 |
| ... | ... | ... | ... | ... |
| 12587 | 417017.587940 | 4.569730e+06 | -100 | 35.940012 |
| 12588 | 417110.834171 | 4.569730e+06 | -100 | 35.901846 |
| 12589 | 417204.080402 | 4.569730e+06 | -100 | 35.360105 |
| 12590 | 416924.341709 | 4.569801e+06 | -100 | 33.214541 |
| 12591 | 417017.587940 | 4.569801e+06 | -100 | 33.599572 |
12592 rows × 4 columns
xyzc_b200max=coordinatesc(DATADIR+'basementknn200confidence.csv',[0,1,2,3])
Now we are ready to draw a figure with the basement surface and the confidence of the prediction of each point where basement is first predicted.
fig=go.Figure()
fig.add_trace(go.Scatter3d(x=xyzc_b200max[0], y=xyzc_b200max[1], z=xyzc_b200max[2],
mode ='markers',
name='knn_basament',
marker = dict(size = 2,
color = xyzc_b200max[3],
colorscale='rainbow',
opacity = 1,
symbol='circle',
showscale=True)
))
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2],
mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
))
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01))
fig.update_layout( title="3D Basement surface Llobregat Delta and confidence index, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Basement_Confidence.html',validate=True, auto_open=False)
'figures/3D_Basement_Confidence.html'
3D figure with the confidence of the basement at each heigh (without the basement surface).
fig=go.Figure()
fig.add_trace(go.Scatter3d(x=x_bas200, y=y_bas200, z=z_bas200,
mode ='markers',
name='knn_basament',
marker = dict(size = 2,
color = conf_bas200,
colorscale='rainbow',
opacity = 1,
symbol='circle',
showscale=True)
))
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2],
mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
#fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
# x=interpolation_knnbasement50cutted[1],
# y=interpolation_knnbasement50cutted[2],
# opacity = 0.5,
# #inensity=interpolation_knnbasementcutted[0],
# colorscale='oryel',
# name='basement surface',
# showscale=False
# ))
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01))
fig.update_layout( title="3D Basement surface Llobregat Delta and confidence index all height, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'Basement_Confidence_Layers.html',validate=True, auto_open=False)
'figures/Basement_Confidence_Layers.html'
Next, we will obtain two figures. The first one shows the gravel predictions colored by confidence. The second one shows the same for the sand predictions.
Confidence list.
conf_gr=np.concatenate([gravels_knn200[i][3] for i in range(21)])
conf_sn=np.concatenate([sands_knn200[i][3] for i in range(21)])
Averages.
gr200_mean=np.mean(conf_gr)
sn200_mean=np.mean(conf_sn)
gr200_mean # Average confidence for gravel predictions
49.33703116232511
sn200_mean # Average confidence for sand predictions
50.18378710802425
x_gr=np.concatenate([gravels_knn200[i][0] for i in range(21)]) # X coordiates of gravel predictions
y_gr=np.concatenate([gravels_knn200[i][1] for i in range(21)]) # Y coordiates of gravel predictions
z_gr=np.concatenate([gravels_knn200[i][2] for i in range(21)]) # Z coordiates of gravel predictions
The following code cell generate the figure with the confidence of the gravel predictions.
fig=go.Figure()
x=x_gr
y=y_gr
z=z_gr
c=conf_gr
fig.add_trace(go.Scatter3d(x=x, y=y, z=z,
mode ='markers',
name='knn_gravels',
marker = dict(size = 2,
color = c,
colorscale='rainbow',
opacity = 1,
symbol='circle',
showscale=True)
))
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2],
mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01))
fig.update_layout( title="3D confidence map for gravels index, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'Gravel_Lithosomes_Confidence_Layers.html',validate=True, auto_open=False)
'figures/Gravel_Lithosomes_Confidence_Layers.html'
Now, we do the same for sand predictions.
x_sn=np.concatenate([sands_knn200[i][0] for i in range(21)])
y_sn=np.concatenate([sands_knn200[i][1] for i in range(21)])
z_sn=np.concatenate([sands_knn200[i][2] for i in range(21)])
fig=go.Figure()
x=x_sn
y=y_sn
z=z_sn
c=conf_sn
fig.add_trace(go.Scatter3d(x=x, y=y, z=z,
mode ='markers',
name='knn_sands',
marker = dict(size = 2,
color = c,
colorscale='rainbow',
opacity = 1,
symbol='circle',
showscale=True)
))
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2],
mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01))
fig.update_layout( title="3D confidence map for sands index, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'Sand_Lithosomes_Confidence_Layers.html',validate=True, auto_open=False)
'figures/Sand_Lithosomes_Confidence_Layers.html'
We shall plot gravel and sand lithosomes like in the previous work https://www.mdpi.com/2077-1312/10/7/986. Then, we will study the confience of the predictions in each lithosome. We shall display tables with the average confidence and draw the lithosomes colored by average confidence.
# knn algorith predictions for a grid of size 50
data_split=[[],[],[],[]]
for h in Heights:
cf=conf(boreholes_data,contour_poly,h,ax,50)
cc=knn(boreholes_data,contour_poly,h,ax,50)
for c in range(4):
x=[cc[0][i,j] for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
y=[cc[1][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
z=[h for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
ccf=[cf[2][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
data_split[c].append([x,y,z,ccf])
gr_knn50=data_split[0] # knn predictions for gravels
sn_knn50=data_split[1] # knn predictions for sands
# We split X, Y, Z and confidence. For gravels
x_gr50=np.concatenate([gr_knn50[i][0] for i in range(21)])
y_gr50=np.concatenate([gr_knn50[i][1] for i in range(21)])
z_gr50=np.concatenate([gr_knn50[i][2] for i in range(21)])
c_gr50=np.concatenate([gr_knn50[i][3] for i in range(21)])
xyzc_gr50=list(zip(x_gr50,y_gr50,z_gr50,c_gr50))
# Now for sands
x_sn50=np.concatenate([sn_knn50[i][0] for i in range(21)])
y_sn50=np.concatenate([sn_knn50[i][1] for i in range(21)])
z_sn50=np.concatenate([sn_knn50[i][2] for i in range(21)])
c_sn50=np.concatenate([sn_knn50[i][3] for i in range(21)])
xyzc_sn50=list(zip(x_sn50,y_sn50,z_sn50,c_sn50))
# We compute the average confidence
gr50_mean=np.mean(c_gr50) # for gravel predictions
sn50_mean=np.mean(c_sn50) # for sand predictions
gr50_mean,sn50_mean
(49.57830307839629, 50.20226163638911)
We compare with those obtained with grid size=200. There is not much difference.
gr200_mean,sn200_mean
(49.33703116232511, 50.18378710802425)
We compute the lithosomes as in https://www.mdpi.com/2077-1312/10/7/986.
# we take the x ,y and z coordinates for the knn predictions
gravels_knn50=[xyc(boreholes_data,contour_poly,ax,0,h,50) for h in Heights]
sands_knn50=[xyc(boreholes_data,contour_poly,ax,1,h,50) for h in Heights]
# Points for gravels.
# p1=(414522,4568893,-115) we are taking max depth = -100, so this lithosome is removed
p2=(418355, 4570004,-80)
p3=(423682,4573012,-60)
p4=(424747,4574333,-65)
p5=(420960, 4571183,-65)
p6=(427019,4578055,-50)
p7=(422096,4577769,-30)
p8=(414901,4570897,-30)
p9=(417173,4571183,-30)
p10=(427777,4579773,-30)
p11=(430427,4579773,-25)
p12=(418309,4573474,-20)
p13=(417931,4582064,-10)
p14=(425126,4577769,0)
gravels_selection1=[#p1,
p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14]
gravels_selection2=[[x[0] for x in gravels_selection1],
[x[1] for x in gravels_selection1],
[x[2] for x in gravels_selection1]]
# Points for sands.
q1=(423611,4570611,-100)
q2=(415658,4568606,-90)
q3=(422096,4570038,-85)
q4=(425883,4572042,-80)
q5=(427398,4576623,-70)
q6=(423990,4575192,-55)
q7=(414144,4568893,-15)
q8=(418688,4570611,-10)
q9=(419445,4574903,-5)
q10=(428155,4576910,-5)
q11=(420203,4572901,0)
q12=(423232,4572329,0)
q13=(422475,4576051,0)
q14=(428155,4574333,0)
q15=(425883,4578628,0)
q16=(427398,4578055,0)
q17=(430806,4580059,0)
sands_selection1=[q1,q2,q3,q4,q5,q6,q7,q8,q9,q10,q11,q12,q13,q14,q15,q16,q17]
sands_selection2=[[x[0] for x in sands_selection1],
[x[1] for x in sands_selection1],
[x[2] for x in sands_selection1]]
# lithosomes for gravels
xygr_knn50=[list(zip(x[0],x[1])) for x in gravels_knn50]
#gravels_lit1=lithosome(p1,xygr_knn50,400,Heights,'grlit1')
#gravels_lit1_b=[x for x in gravels_lit1[0] if x[2]< -35]
#h1_b = ConvexHull(gravels_lit1_b)
#gravels_lit1b=[h1_b.points,h1_b.vertices,h1_b.simplices,'grlit1']
gravels_lit2=lithosome(p2,xygr_knn50,400,Heights,'grlit1')
gravels_lit3=lithosome(p3,xygr_knn50,400,Heights,'grlit2')
gravels_lit4=lithosome(p4,xygr_knn50,400,Heights,'grlit3')
gravels_lit5=lithosome(p5,xygr_knn50,400,Heights,'grlit4')
gravels_lit6=lithosome(p6,xygr_knn50,400,Heights,'grlit5')
gravels_lit7=lithosome(p7,xygr_knn50,400,Heights,'grlit6')
gravels_lit8=lithosome(p8,xygr_knn50,400,Heights,'grlit7')
gravels_lit9=lithosome(p9,xygr_knn50,400,Heights,'grlit8')
gravels_lit10=lithosome(p10,xygr_knn50,400,Heights,'grlit9')
gravels_lit11=lithosome(p11,xygr_knn50,500,Heights,'grlit10')
gravels_lit12=lithosome(p12,xygr_knn50,500,Heights,'grlit11')
gravels_lit13=lithosome(p13,xygr_knn50,400,Heights,'grlit12')
gravels_lit14=lithosome(p14,xygr_knn50,400,Heights,'grlit13')
# List of gravels lithosomes
list_grlithosomes=[#gravels_lit1b,
gravels_lit2,gravels_lit3,gravels_lit4,
gravels_lit5,gravels_lit6,gravels_lit7,gravels_lit8,
gravels_lit9,gravels_lit10,gravels_lit11,gravels_lit12,
gravels_lit13,gravels_lit14]
# lithosome for sands
xysnd_knn=[list(zip(x[0],x[1])) for x in sands_knn50]
snd_lit1=lithosome(q1,xysnd_knn,400,Heights,'sndlit1')
snd_lit1_b=[x for x in snd_lit1[0] if x[2]< -89]
ha1_1b = ConvexHull(snd_lit1_b)
snd_lit1b=[ha1_1b.points,ha1_1b.vertices,ha1_1b.simplices,'sndlit1']
snd_lit2=lithosome(q2,xysnd_knn,400,Heights,'sndlit2')
snd_lit2_b=[x for x in snd_lit2[0] if x[2]<-70]
ha2_1b = ConvexHull(snd_lit2_b)
snd_lit2b=[ha2_1b.points,ha2_1b.vertices,ha2_1b.simplices,'sndlit2']
snd_lit3=lithosome(q3,xysnd_knn,400,Heights,'sndlit3')
snd_lit3_b=[x for x in snd_lit3[0] if x[2]< -70 and -90<x[2]]+[
[ 4.20960571e+05, 4.56975163e+06, -9.00000000e+01],
[ 4.21339265e+05, 4.56946531e+06, -9.00000000e+01],
[ 4.21339265e+05, 4.56975163e+06, -9.00000000e+01],
[ 4.21339265e+05, 4.57003796e+06, -9.00000000e+01],
[ 4.21717959e+05, 4.56975163e+06, -9.00000000e+01],
[ 4.21717959e+05, 4.57003796e+06, -9.00000000e+01],
[ 4.21717959e+05, 4.57032429e+06, -9.00000000e+01],
[ 4.21717959e+05, 4.57061061e+06, -9.00000000e+01],
[ 4.22096653e+05, 4.56975163e+06, -9.00000000e+01],
[ 4.22096653e+05, 4.57003796e+06, -9.00000000e+01],
[ 4.22096653e+05, 4.57032429e+06, -9.00000000e+01],
[ 4.22096653e+05, 4.57061061e+06, -9.00000000e+01],
[ 4.22475347e+05, 4.57003796e+06, -9.00000000e+01],
[ 4.22475347e+05, 4.57032429e+06, -9.00000000e+01],
[ 4.22475347e+05, 4.57061061e+06, -9.00000000e+01]]
ha3_1b = ConvexHull(snd_lit3_b)
snd_lit3b=[ha3_1b.points,ha3_1b.vertices,ha3_1b.simplices,'sndlit3']
snd_lit4=lithosome(q4,xysnd_knn,400,Heights,'sndlit4')
snd_lit4_b=[x for x in snd_lit4[0] if x[2]< -10 and x[2]>-90]+[
[ 4.25504898e+05, 4.57118327e+06, -9.00000000e+01],
[ 4.25504898e+05, 4.57146959e+06, -9.00000000e+01],
[ 4.25504898e+05, 4.57175592e+06, -9.00000000e+01],
[ 4.25504898e+05, 4.57204224e+06, -9.00000000e+01],
[ 4.25504898e+05, 4.57232857e+06, -9.00000000e+01],
[ 4.25504898e+05, 4.57261490e+06, -9.00000000e+01],
[ 4.25883592e+05, 4.57232857e+06, -9.00000000e+01],
[ 4.25883592e+05, 4.57261490e+06, -9.00000000e+01],
[ 4.25883592e+05, 4.57290122e+06, -9.00000000e+01],
[ 4.25883592e+05, 4.57318755e+06, -9.00000000e+01],
[ 4.25883592e+05, 4.57347388e+06, -9.00000000e+01],
[ 4.26262286e+05, 4.57175592e+06, -9.00000000e+01],
[ 4.26262286e+05, 4.57204224e+06, -9.00000000e+01],
[ 4.26262286e+05, 4.57232857e+06, -9.00000000e+01],
[ 4.26262286e+05, 4.57261490e+06, -9.00000000e+01],
[ 4.26262286e+05, 4.57290122e+06, -9.00000000e+01],
[ 4.26262286e+05, 4.57318755e+06, -9.00000000e+01],
[ 4.2664098e+05, 4.5726149e+06, -9.0000000e+01],
[ 4.26640980e+05, 4.57290122e+06, -9.00000000e+01]]
ha4_1b = ConvexHull(snd_lit4_b)
snd_lit4b=[ha4_1b.points,ha4_1b.vertices,ha4_1b.simplices,'sndlit4']
snd_lit5=lithosome(q5,xysnd_knn,400,Heights,'sndlit5')
snd_lit5_b=[x for x in snd_lit5[0] if x[2]<-50]
ha5_1b = ConvexHull(snd_lit5_b)
snd_lit5b=[ha5_1b.points,ha5_1b.vertices,ha5_1b.simplices,'sndlit5']
snd_lit6=lithosome(q6,xysnd_knn,400,Heights,'sndlit6')
snd_lit6_b=[x for x in snd_lit6[0] if x[2]<-50]
ha6_1b = ConvexHull(snd_lit6_b)
snd_lit6b=[ha6_1b.points,ha6_1b.vertices,ha6_1b.simplices,'sndlit6']
snd_lit7=lithosome(q7,xysnd_knn,400,Heights,'sndlit7')
snd_lit7_b=[x for x in snd_lit7[0] if x[2]>-30]
ha7_1b = ConvexHull(snd_lit7_b)
snd_lit7b=[ha7_1b.points,ha7_1b.vertices,ha7_1b.simplices,'sndlit8']
snd_lit8=lithosome(q8,xysnd_knn,400,Heights,'sndlit8')
snd_lit8_b=[x for x in snd_lit8[0] if x[2]>-30 and x[2]<-4]
ha8_1b = ConvexHull(snd_lit8_b)
snd_lit8b=[ha8_1b.points,ha8_1b.vertices,ha8_1b.simplices,'sndlit8']
snd_lit9=lithosome(q9,xysnd_knn,400,Heights,'sndlit9')
snd_lit9_b=[x for x in snd_lit9[0] if x[2]>-25]
ha9_1b = ConvexHull(snd_lit9_b)
snd_lit9b=[ha9_1b.points,ha9_1b.vertices,ha9_1b.simplices,'sndlit9']
snd_lit10=lithosome(q10,xysnd_knn,400,Heights,'sndlit10')
snd_lit10_b=[x for x in snd_lit10[0] if x[2]>-35]
ha10_2b = ConvexHull(snd_lit10_b)
snd_lit10b=[ha10_2b.points,ha10_2b.vertices,ha10_2b.simplices,'sndlit10']
snd_lit11=lithosome(q11,xysnd_knn,400,Heights,'sndlit11')
snd_lit11_b=[x for x in snd_lit11[0] if x[2]>-36]
ha11_1b = ConvexHull(snd_lit11_b)
snd_lit11b=[ha11_1b.points,ha11_1b.vertices,ha11_1b.simplices,'sndlit11']
snd_lit12=lithosome(q12,xysnd_knn,400,Heights,'sndlit12')
snd_lit12_b=[x for x in snd_lit12[0] if x[2]>-37 and x[1]<4576053]
ha12_2b = ConvexHull(snd_lit12_b)
snd_lit12b=[ha12_2b.points,ha12_2b.vertices,ha12_2b.simplices,'sndlit12']
snd_lit13=lithosome(q13,xysnd_knn,400,Heights,'sndlit13')
snd_lit13_b=[x for x in snd_lit13[0] if x[2]>-27]
ha13_3b = ConvexHull(snd_lit13_b)
snd_lit13b=[ha13_3b.points,ha13_3b.vertices,ha13_3b.simplices,'sndlit13']
snd_lit14=lithosome(q14,xysnd_knn,400,Heights,'sndlit14')
snd_lit14_b=[x for x in snd_lit14[0] if x[2]>-57]
ha14_4b = ConvexHull(snd_lit14_b)
snd_lit14b=[ha14_4b.points,ha14_4b.vertices,ha14_4b.simplices,'sndlit14']
snd_lit15=lithosome(q15,xysnd_knn,400,Heights,'sndlit15')
snd_lit15_b=[x for x in snd_lit15[0] if x[1]>4578055 ]
ha15_5b = ConvexHull(snd_lit15_b)
snd_lit15b=[ha15_5b.points,ha15_5b.vertices,ha15_5b.simplices,'sndlit15']
snd_lit16=lithosome(q16,xysnd_knn,400,Heights,'sndlit16')
snd_lit16_b=[x for x in snd_lit16[0] if x[2]>-60 and x[1]>4576020 and x[1]<4578100 and x[0]>425883]
ha16_6b = ConvexHull(snd_lit16_b)
snd_lit16b=[ha16_6b.points,ha16_6b.vertices,ha16_6b.simplices,'sndlit16']
snd_lit17=lithosome(q17,xysnd_knn,500,Heights,'sndlit17')
# List of sands lithosomes
list_sndlithosomesb=[snd_lit1b,snd_lit2b,snd_lit3b,snd_lit4b,snd_lit5b,
snd_lit6b,snd_lit7b,snd_lit8b,snd_lit9b,snd_lit10b,
snd_lit11b,snd_lit12b,snd_lit13b,snd_lit14b,snd_lit15b,
snd_lit16b,snd_lit17]
list_sndlithosomes=[snd_lit1,snd_lit2,snd_lit3,snd_lit4,snd_lit5,
snd_lit6,snd_lit7,snd_lit8,snd_lit9,snd_lit10,
snd_lit11,snd_lit12,snd_lit13,snd_lit14,snd_lit15,
snd_lit16,snd_lit17]
First, we reshape the data.
len(list_grlithosomes[0])
5
len(list_sndlithosomes[0])
5
zip_list_grlithosomes=[list(zip(x[4][0],x[4][1],x[4][2])) for x in list_grlithosomes]
We group the gravel predictions by lithosome.
# data points
list_knn_in_grlit=[[x for x in xyzc_gr50 if ((x[0],x[1],x[2]) in y)] for y in zip_list_grlithosomes]
# confidence values
list_of_c_grlit=[[x[3] for x in y] for y in list_knn_in_grlit]
# heights
z_grlit=[[x[2] for x in y] for y in list_knn_in_grlit]
We work out the average confidence and maximum and minimum heights for each lithosome.
The auxiliar function mean_c2 computes the means.
def mean_c2(lit):
zlit=list(zip(lit[4][0],lit[4][1],lit[4][2]))
knlit=[x for x in xyzc_gr50 if ((x[0],x[1],x[2]) in zlit)]
list_of_c_grlit=[x[3] for x in knlit]
return np.mean(list_of_c_grlit)
m2=[mean_c2(x) for x in list_grlithosomes]
The following lines perform the same computation as the function mean_c2.
means_cgrlit=[np.mean(x) for x in list_of_c_grlit] # For confidence value
means_zgrlit=[np.mean(x) for x in z_grlit] # For Z coordinate
We check that the same result is obtained.
m2==means_cgrlit
True
We compute the maximum and minimum height for each lithosome.
max_zgrlit=[max(x) for x in z_grlit]
min_zgrlit=[min(x) for x in z_grlit]
We create the table.
lit_gr_names=['grlit'+str(i+1) for i in range(13)]
d1={'graves lithosomes':lit_gr_names,
'confidence in mean':means_cgrlit,
'average height':means_zgrlit,
'max height':max_zgrlit,
'min height':min_zgrlit
}
df = pd.DataFrame(d1)
df.to_csv(DATADIR+'gr_lit_table.csv', index=False, header=True)
df
| graves lithosomes | confidence in mean | average height | max height | min height | |
|---|---|---|---|---|---|
| 0 | grlit1 | 50.174238 | -72.397590 | -40 | -100 |
| 1 | grlit2 | 51.354275 | -51.391417 | -40 | -60 |
| 2 | grlit3 | 51.747555 | -55.541150 | -45 | -70 |
| 3 | grlit4 | 49.055077 | -56.647059 | -50 | -65 |
| 4 | grlit5 | 50.939566 | -45.714286 | -35 | -50 |
| 5 | grlit6 | 49.723211 | -28.293173 | 0 | -40 |
| 6 | grlit7 | 47.790273 | -29.807692 | -10 | -40 |
| 7 | grlit8 | 43.996622 | -29.513889 | -25 | -35 |
| 8 | grlit9 | 53.124705 | -19.677419 | -10 | -25 |
| 9 | grlit10 | 28.545459 | -20.000000 | -15 | -25 |
| 10 | grlit11 | 45.200411 | -17.500000 | -15 | -20 |
| 11 | grlit12 | 49.076593 | -12.685590 | 0 | -25 |
| 12 | grlit13 | 50.428450 | -25.000000 | 0 | -40 |
We can compare the average confidence when taking a grid of size 50 or 200.
gr50_mean # Average confidence for gravel predictions with grid size 50
49.57830307839629
gr200_mean # Average confidence for gravel predictions with grid size 200
49.33703116232511
len(list_sndlithosomes)
17
zip_list_snlithosomes=[list(zip(x[4][0],x[4][1],x[4][2])) for x in list_sndlithosomes]
# data points
list_knn_in_snlit=[[x for x in xyzc_sn50 if ((x[0],x[1],x[2]) in y)] for y in zip_list_snlithosomes]
# We group the sand predictions by lithosome
# Confidence value
list_of_c_snlit=[[x[3] for x in y] for y in list_knn_in_snlit]
# Heights
z_snlit=[[x[2] for x in y] for y in list_knn_in_snlit]
# We work out the average confidence and maximum and minimum heights for each lithosome
means_csnlit=[np.mean(x) for x in list_of_c_snlit]
means_zsnlit=[np.mean(x) for x in z_snlit]
max_zsnlit=[max(x) for x in z_snlit]
min_zsnlit=[min(x) for x in z_snlit]
# Table
lit_sn_names=['snlit'+str(i+1) for i in range(17)]
d2={'sand lithosomes':lit_sn_names,
'confidence in mean':means_csnlit,
'cota mean':means_zsnlit,
'cota max':max_zsnlit,
'cota min':min_zsnlit
}
dff = pd.DataFrame(d2)
dff.to_csv(DATADIR+'sn_lit_table.csv', index=False, header=True)
dff
| sand lithosomes | confidence in mean | cota mean | cota max | cota min | |
|---|---|---|---|---|---|
| 0 | snlit1 | 50.516362 | -83.949275 | 0 | -100 |
| 1 | snlit2 | 54.352117 | -35.000000 | 0 | -90 |
| 2 | snlit3 | 53.833335 | -25.512465 | 0 | -90 |
| 3 | snlit4 | 52.791327 | -39.709507 | 0 | -100 |
| 4 | snlit5 | 41.704941 | -65.909091 | 0 | -80 |
| 5 | snlit6 | 53.765592 | -12.724057 | 0 | -60 |
| 6 | snlit7 | 54.228871 | -24.390244 | 0 | -90 |
| 7 | snlit8 | 53.585398 | -8.754941 | 0 | -40 |
| 8 | snlit9 | 50.545163 | -10.555556 | -5 | -20 |
| 9 | snlit10 | 46.070302 | -49.139785 | -5 | -80 |
| 10 | snlit11 | 45.906706 | -17.142857 | 0 | -40 |
| 11 | snlit12 | 54.897272 | -10.053763 | -5 | -55 |
| 12 | snlit13 | 42.450885 | -40.594059 | 0 | -60 |
| 13 | snlit14 | 52.951239 | -16.250000 | 0 | -60 |
| 14 | snlit15 | 54.124507 | -4.559322 | 0 | -10 |
| 15 | snlit16 | 53.311568 | -9.187117 | 0 | -60 |
| 16 | snlit17 | 29.756034 | -6.562500 | 0 | -15 |
sn50_mean # Average confidence for sand predictions with grid size 50
50.20226163638911
sn200_mean # Average confidence for sand predictions with grid size 200
50.18378710802425
# Lithosome data for gravels
data_grlithosomes=[data_lithosome(x[0],x[2],x[3],0,1,'skyblue') for x in list_grlithosomes]
# Lithosome data for sands
data_sndlithosomes=[data_lithosome(x[0],x[2],x[3],0,1,'lemonchiffon') for x in list_sndlithosomesb]
# lithosomes for gravels
dat=data_grlithosomes
fig=go.Figure(data=dat)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout( title="3D graves lithosomes Llobregat Delta, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes.html',validate=True, auto_open=False)
'figures/3D_grlithosomes.html'
# lithosomes for sands
dat=data_sndlithosomes
fig=go.Figure(data=dat)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout( title="3D sands lithosomes Llobregat Delta, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes.html',validate=True, auto_open=False)
'figures/3D_grlithosomes.html'
# lithosomes for gravels and sands
dat=data_grlithosomes+data_sndlithosomes
fig=go.Figure(data=dat)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout( title="3D graves and sands lithosomes Llobregat Delta, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_lithosomes.html',validate=True, auto_open=False)
'figures/3D_lithosomes.html'
With basement.
# lithosomes for gravels and sands with basement
dat=data_grlithosomes+data_sndlithosomes
fig=go.Figure(data=dat)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
))
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout( title="3D graves and sands lithosomes Llobregat Delta with basement surface, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Lithosomes_And_Basement_LRD.html',validate=True, auto_open=False)
'figures/3D_Lithosomes_And_Basement_LRD.html'
We modify the function data_lithosome to employ different colors.
def data_lithosome2(vhull,shull,name,alpha,opc,col):
return go.Mesh3d(x=vhull[:, 0],y=vhull[:, 1], z=vhull[:, 2],
name=name,
showlegend=True,
colorbar_title=name,
color=col,
colorscale='rainbow',
i=shull[:, 0],
j=shull[:, 1],
k=shull[:, 2],
opacity=opc,
alphahull=alpha,
showscale=True
)
# Data for gravels
data_grlithosomes_conf=[data_lithosome2(x[0],x[2],x[3],0,1,mean_c2(x)) for x in list_grlithosomes]
# the complete gravels lithosome data
fig=go.Figure(data=data_grlithosomes_conf)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
))
fig.update_layout( title="3D graves lithosomes Llobregat Delta colored by confidence, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes_con.html',validate=True, auto_open=False)
'figures/3D_grlithosomes_con.html'
len(list_sndlithosomesb)
17
# Data for sands
data_snlithosomes_conf=[data_lithosome2(list_sndlithosomesb[i][0],
list_sndlithosomesb[i][2],
list_sndlithosomesb[i][3],
0,1,mean_c2(list_sndlithosomes[i])) for i in range(17)]
fig=go.Figure(data=data_snlithosomes_conf)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
))
fig.update_layout( title="3D sand lithosomes Llobregat Delta colored by confidence, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_snlithosomes_con.html',validate=True, auto_open=False)
'figures/3D_snlithosomes_con.html'
Together with basement.
# Data for gravels with opacity 0.3
data_grlithosomes_conf05=[data_lithosome2(x[0],x[2],x[3],0,0.5,mean_c2(x)) for x in list_grlithosomes]
dat=data=data_snlithosomes_conf+data_grlithosomes_conf
fig=go.Figure(data=dat)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
))
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)))
fig.update_layout( title="3D gravels and sand lithosomes Llobregat Delta colored by confidence, Z scale is x 50.",
scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_lithosomes_con.html',validate=True, auto_open=False)
'figures/3D_lithosomes_con.html'
We will now display two figures side by side. One with the lithosomes in a single color and another with lithosomes colored by the average confidence.
We need to import subplots from plotly.
from plotly.subplots import make_subplots
Without basement.
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
subplot_titles=('graves lithosomes','graves lithosomes colored by confidence')
)
for i in range(len(data_grlithosomes)):
fig.add_trace(data_grlithosomes[i],row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
legendgroup = '1',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=1)
#fig.update_layout(legend=dict(x=0,y=1))
for i in range(len(data_grlithosomes_conf)):
fig.add_trace(data_grlithosomes_conf[i],row=1, col=2)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=2)
fig.update_traces(showlegend=True)
fig.update_layout( title="3D graves lithosomes comparation Llobregat Delta, Z scale is x 50.",
scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
legend=dict(x=1,y=0),
height=600, showlegend=True)
#fig.update_layout(scene2_aspectmode='manual',
# scene2_aspectratio=dict(x=1, y=1, z=0.5),
# legend=dict(x=1,y=0)
# )
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes_comparation_nobasement.html',validate=True, auto_open=False)
'figures/3D_grlithosomes_comparation_nobasement.html'
With basement.
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
subplot_titles=('graves lithosomes','graves lithosomes colored by confidence')
)
for i in range(len(data_grlithosomes)):
fig.add_trace(data_grlithosomes[i],row=1, col=1)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
), row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
legendgroup = '1',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=1)
#fig.update_layout(legend=dict(x=0,y=1))
for i in range(len(data_grlithosomes_conf)):
fig.add_trace(data_grlithosomes_conf[i],row=1, col=2)
#fig.add_trace(marks_data[0],row=1, col=2)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
), row=1, col=2)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=2)
fig.update_traces(showlegend=True)
fig.update_layout( title="3D graves lithosomes comparation Llobregat Delta, Z scale is x 50.",
scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
legend=dict(x=1,y=0),
height=600, showlegend=True)
#fig.update_layout(scene2_aspectmode='manual',
# scene2_aspectratio=dict(x=1, y=1, z=0.5),
# legend=dict(x=1,y=0)
# )
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Gravel_Lithosomes_Confidence.html',validate=True, auto_open=False)
'figures/3D_Gravel_Lithosomes_Confidence.html'
Without basement.
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
subplot_titles=('sands lithosomes','sands lithosomes colored by confidence')
)
for i in range(len(data_sndlithosomes)):
fig.add_trace(data_sndlithosomes[i],row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=1)
for i in range(len(data_snlithosomes_conf)):
fig.add_trace(data_snlithosomes_conf[i],row=1, col=2)
#fig.add_trace(marks_data[0],row=1, col=2)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=2)
#fig.update_layout(legend=dict(
# yanchor="top",
# y=0.99,
# xanchor="right",
#x=0.01))
fig.update_layout( title="3D sands lithosomes comparation Llobregat Delta, Z scale is x 50.",
scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
height=600, showlegend=True)
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_snlithosomes_comparation_nobasement.html',validate=True, auto_open=False)
'figures/3D_snlithosomes_comparation_nobasement.html'
With basement.
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
subplot_titles=('sands lithosomes','sands lithosomes colored by confidence')
)
for i in range(len(data_sndlithosomes)):
fig.add_trace(data_sndlithosomes[i],row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
), row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=1)
for i in range(len(data_snlithosomes_conf)):
fig.add_trace(data_snlithosomes_conf[i],row=1, col=2)
#fig.add_trace(marks_data[0],row=1, col=2)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
x=interpolation_knnbasement50cutted[1],
y=interpolation_knnbasement50cutted[2],
opacity = 0.5,
#inensity=interpolation_knnbasementcutted[0],
colorscale='oryel',
name='basement surface',
showscale=False
), row=1, col=2)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
line_width=5,
name='Delta Contour',
marker = dict(
size = 4,
color = 'black'
)), row=1, col=2)
#fig.update_layout(legend=dict(
# yanchor="top",
# y=0.99,
# xanchor="right",
#x=0.01))
fig.update_layout( title="3D sands lithosomes comparation Llobregat Delta, Z scale is x 50.",
scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
#aspectmode='auto',
#aspectratio=dict(x=2, y=2, z=1)
),
height=600, showlegend=True)
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Sand_Lithosomes_Confidence.html',validate=True, auto_open=False)
'figures/3D_Sand_Lithosomes_Confidence.html'